#include <petscksp.h>
#include <assert.h>

int main(int argc,
	 char** argv)
{
    Vec            x, b, u;      /* approx solution, RHS, exact solution */
    Mat            matA;            /* linear system matrix */
    KSP            ksp;         /* linear solver context */
    PC             pc;           /* preconditioner context */
    PetscErrorCode ierr;

    /**
     * Reads CRS format into PETSC SEQAIJ format
     */
    PetscInt nr;
    PetscInt nv;

    /* read case0.dat to matrix dimensions */
    FILE* finCaseDat0;
    finCaseDat0 = fopen("case.dat0", "r");
    fscanf(finCaseDat0, "%d", &nr);
    fscanf(finCaseDat0, "%d", &nv);
    assert(nr == 9);
    assert(nv == 5);
    close(finCaseDat0);
    PetscInt n = nv*nr;
    
    /* allocat matA */
    // PetscReal* matA = calloc(n*n, sizeof(PetscReal));

    ierr = MatCreate(PETSC_COMM_WORLD,
		     &matA);
    CHKERRQ(ierr);
    ierr = MatSetSizes(matA,
		       PETSC_DECIDE,
		       PETSC_DECIDE,
		       n,
		       n);
    CHKERRQ(ierr);
    ierr = MatSetFromOptions(matA);
    CHKERRQ(ierr);
    
    FILE* finV;
    FILE* finR;
    FILE* finC;

    finV = fopen("matrix_coef.dat","r");
    finR = fopen("row_ptr.dat", "r");
    finC = fopen("col_ind.dat", "r");

    PetscInt row = 0;   /* row counter */
    PetscInt rs;  /* start index in row_ptr.dat */
    PetscInt re;  /* end index in row_ptr.dat */
    fscanf(finR, "%d", &rs);
    fscanf(finR, "%d", &re);
    while(feof(finR))
    {
	PetscInt range = re - rs;
	PetscInt i;
	for (i = 0; i < range; ++i)
	{
	    assert(!feof(finC));
	    PetscInt col;
	    fscanf(finC, "%d", &col);
	    
	    PetscInt j;
	    for (j = 0; j < nv; ++j)
	    {
		PetscInt cols[nv];
		PetscScalar vals[nv];
		PetscInt k;
		for (k = 0; k < nv; k++)
		{
		    assert(!feof(finV));
		    PetscReal v;
		    fscanf(finV, "%le", &v);
		    cols[k] = col*nv+k;
		    vals[k] = v;
		}
		PetscInt itmp = row*nv+j;
		ierr = MatSetValues(matA, 1, &itmp, 5, cols, vals, INSERT_VALUES); 
	    }
	}
	
	/* update to another row */
	rs = re;
	assert(!feof(finR));
	re = fscanf(finR, "%d", &re);
	++row;
    }

    ierr = MatView(matA, PETSC_VIEWER_STDOUT_SELF);
    CHKERRQ(ierr);
    
    PetscPrintf(PETSC_COMM_WORLD, "%s", "...done\n");
    CHKERRQ(ierr);

    return 0;
}
